#ifndef EVECTOR_HPP_
#define EVECTOR_HPP_

/* The following code example is taken from the book
 * "C++ Templates - The Complete Guide"
 * by David Vandevoorde and Nicolai M. Josuttis, Addison-Wesley, 2002
 *
 * (C) Copyright David Vandevoorde and Nicolai M. Josuttis 2002.
 * Permission to copy, use, modify, sell and distribute this software
 * is granted provided this copyright notice appears in all copies.
 * This software is provided "as is" without express or implied
 * warranty, and with no claim as to its suitability for any purpose.
 */
/* modified by ofgr 21.03.2008 */
#include <stddef.h>
#include <iostream>
#include <cassert>
#include "svector.hpp"
#include "usertypes.hpp"


// *****************************************
// class for objects that represent scalars
// *****************************************
template <typename T>
class A_Scalar {
  private:
    T const& s;  // value of the scalar

  public:
    // constructor initializes value
    A_Scalar (T const& v)
     : s(v) {
    }

    // for index operations the scalar is the value of each element
    T operator[] (myindex) const {
        return s;
    }

    // scalars have zero as size
    myindex size() const {
        return 0;
    };
};

// *****************************************
// class for objects that represent vectors
// *****************************************
template <typename T, typename Rep = SVector<T> >
class EVector 
{
private:
  Rep expr_rep;   // (access to) the data of the vector
  
public:
	  double GetNorm2();
	    T ScalarProd(EVector<T,Rep> other);
	// create vector with initial size

  explicit EVector (myindex s)
    : expr_rep(s) {}

  // create vector from possible representation
  EVector (Rep const& rb)
    : expr_rep(rb) {}

  // assignment operator for same type
  EVector& operator= (EVector const& b) 
  { 
    assert(size()==b.size());
    for (myindex idx = 0; idx<b.size(); ++idx)
      expr_rep[idx] = b[idx];
    
    return *this;
  }

  // assignment operator for vectors of different type
  template<typename T2, typename Rep2>
  EVector& operator= (EVector<T2, Rep2> const& b) 
  { 
    assert(size()==b.size());
    for (myindex idx = 0; idx<b.size(); ++idx)
      expr_rep[idx] = b[idx];
    
    return *this;
  }

  // set vector zero
  void setZero ()
  {
    for (myindex idx = 0; idx<size(); ++idx)
      expr_rep[idx] = 0;
  }

  // += operator for vectors of different type
  template<typename T2, typename Rep2>
  EVector& operator+= (EVector<T2, Rep2> const& b) 
  { 
    assert(size()==b.size());
    for (myindex idx = 0; idx<b.size(); ++idx)
      expr_rep[idx] += b[idx];
    
    return *this;
  }

  // += operator for vector and scalar of the same type
  EVector& operator+= (A_Scalar<T> const& b) 
  { 
    for (myindex idx = 0; idx<size(); ++idx)
      expr_rep[idx] += b[idx];
    
    return *this;
  }

  // -= operator for vectors of different type
  template<typename T2, typename Rep2>
  EVector& operator-= (EVector<T2, Rep2> const& b) 
  { 
    assert(size()==b.size());
    for (myindex idx = 0; idx<b.size(); ++idx)
      expr_rep[idx] -= b[idx];
    
    return *this;
  }

 // -= operator for vector and scalar of the same type
  EVector& operator-= (A_Scalar<T> const& b) 
  { 
    for (myindex idx = 0; idx<size(); ++idx)
      expr_rep[idx] -= b[idx];
    
    return *this;
  }

  // *= operator for vectors of different type
  template<typename T2, typename Rep2>
  EVector& operator*= (EVector<T2, Rep2> const& b) 
  { 
    assert(size()==b.size());
    for (myindex idx = 0; idx<b.size(); ++idx)
      expr_rep[idx] *= b[idx];
    
    return *this;
  }

  // *= operator for vector and scalar of the same type
  EVector& operator*= (A_Scalar<T> const& b) 
  { 
    for (myindex idx = 0; idx<size(); ++idx)
      expr_rep[idx] *= b[idx];
    
    return *this;
  }  



  // size is size of represented data
  myindex size() const 
  {
    return expr_rep.size();
  }

  // index operator for constants and variables
  T operator[] (myindex idx) const 
  {
    //assert(idx<size());
    return expr_rep[idx];
  }
  
  T& operator[] (myindex idx) 
  {
    //assert(idx<size());
    return expr_rep[idx];
  }

  // return what the vector currently represents
  Rep const& rep() const 
  {
    return expr_rep;
  }
  Rep& rep() 
  {
    return expr_rep;
  }
  // print 
  void print () 
  { 
    for (myindex idx = 0; idx<size(); ++idx)
      std::cout << expr_rep[idx] << std::endl;
  }


};
 template<typename T2, typename Rep2> double EVector<T2,Rep2>::GetNorm2()
{
    double sum = 0;
    for(int i = 0; i < size(); i++)
        sum +=expr_rep[i] *expr_rep[i];
    return sqrt(sum);
}
template<typename T>
std::ostream& operator<< (std::ostream& s, const EVector<T,SVector<T> >& b)
{
  for (myindex i=0; i<b.size(); i++)
    s << " " << b[i];

  return (s);
}
template<typename T2, typename Rep2> T2 EVector<T2,Rep2>::ScalarProd(EVector<T2,Rep2> other)
{
    if(size() != other.size())
    {
        cout << "Dimensions must agree!";
        throw;
    }
    else
    {
        T2 sum = 0;
        for(int i = 0; i < size(); i++)
            sum += expr_rep[i] * other.expr_rep[i];
        return sum;
    }
}
#endif // EVECTOR_
